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ABS^ACT 

\j/ 

v  Our  study  of  the  interaction  of  millimeter  wave  fields  with  biological  systems 
has  concentrated  on  Davydov's  mw^el  of  sol  1 ton  formation  In  alpha-helices.  A 
qualitative  understanding  has  been  obtained  for  the  oscillatory  modes  of  proteins  In 
the  millimeter  and  submllllmeter  regions  of  the  electromagnetic  spectrum. 
Specifically: 

o  ''i  Four  separate  groups  of  frequencies  are  found  to  exist  in  short  alpha- 
helices,  ranging  from  200  gigahertz  to  6  terahertz. 

o  In  long  alpha-hel 1ces»  sollton  trapping  Is  predicted  for  multiquanta 
excitation  at  zero  temperature*  while  the  existence  of  solltons  at  room 
temperature  Is  called  into  question. 

o  A  molecular  orbital  calculation  of  the  excl ton-phonon  coupling  constant  In 
the  formamlde  dimer  suggests  that  the  value  of  this  critical  constant 
still  needs  to  be  determined. 

The  results  point  up  the  need  for  further  molecular  orbital  calculations  and 
experimental  verification  In  order  to  understand  the  interaction  of  millimeter 
waves  with  biological  systems.  ^ 
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FINAL  REPORT 


A  STUDY  OF  THE  INTERACTION  OF  MILLIMETER  WAVE  FIELDS  WITH  BIOLOGICAL  SYSTEMS 


I.  INTRODUCTION 

Recent  experimental  and  theoretical  work  Indicates  that  some  biological  systems 
exhibit  sensitivities  to  electromagnetic  radiation  In  the  millimeter  and 
submil llmeter  wavelengths.  These  fields,  even  at  very  low  Intensities,  have  been 
reported  to  produce  profound  and  marked  alterations  In  a  variety  of  organisms. 
Specific  changes  Include  enhancement  and  suppression  of  growth  rates  In  Saccaromvces 
cerevlslae  (1,2),  modification  of  genetic  materials  In  Drosophila,  (3),  and 
retardation  of  free  radical  reaction  rates  In  Rhodopseudomonas  Sohaeroldes  (4). 

In  Industrial  or  military  environments  a  significant  number  of  personnel  are 
exposed  to  electromagnetic  fields  comprised  of  two  or  more  frequencies.  The  need  to 
account  for  the  effects  of  compound  fields  places  a  special  burden  on  biological 
theory.  Research  with  simple  sinusoidal  fields  Indicates  that  the  mechanisms  are 
nonlinear,  that  both  amplitude  and  frequency  are  significant  In  accounting  for 
biological  effects.  Granting  that  both  amplitude  and  frequency  windows  are  present 
for  simple  sinusoidal  fields,  the  possibilities  for  compound  windows  In 
superpositions  of  two  or  more  sinusoidal  fields  Increases  the  number  of  possible 
Interactions  beyond  what  Is  reasonable  for  case-by-case  testing.  Responses  to  the 
extremely  low  frequency  (ELF)  range  (0.1  to  20000  Hz)  and  the  microwave  to 
submllllmeter  range  (300  MHz  to  100  GHz)  have  particular  interest  because  of  the 
challenge  to  current  theories  of  Interactions  between  electromagnetic  radiation  and 
living  tissue. 

It  Is  necessary  to  address  the  Issues  from  the  perspective  of  the  physical 
sciences  In  order  to  establish  which  biological  processes  are  relevant  to  the  EMF 
sensitivities  of  living  systems.  Adey  (5)  has  recently  proposed  that  the  cell 
membrane  Is  the  most  likely  site  of  Interaction  for  the  electromagnetic  field. 

Adey's  model  of  transmembrane  coupling  between  extracellular  fields  and 
Intracellular  responses  leaves  unspecified  several  requirements  for  an  amplifying 
mechanism  between  Initial  cell  surface  events  and  resulting  enzymatic  responses. 

The  first  requirement  Is  for  an  amplification  of  the  weak  extracellular 
electromagnetic  fields  sufficiently  to  perturb  binding  events  at  the  cell  surface. 
The  second  Is  for  a  means  to  pass  a  signal  through  the  barrier  presented  by  the  cell 
membrane.  The  third  Is  for  a  source  of  energy  to  make  up  the  difference  between 
molecular  binding  energies  available  at  the  cell  surface  and  energies  released  by 
subsequent  Intracellular  events. 

During  the  past  year  we  have  Investigated  theoretical  mechanisms  for  energy 
transfer  between  ambient  fields  and  Ions  or  molecules  at  the  cell  surface  to 
proteins  in  the  Interior  of  the  cell.  Although  current  research  Is  proceeding 
toward  an  unambiguous  Identification  of  cell  membrane  systems  Involved  In  the 
transduction  of  weak  electrical  fields  (see  Adey,  (6))  data  on  the  specific 
proteins  Involved  Is  presently  unavailable.  Knowledge  of  protein  structure  Is 
essential  to  a  study  of  protein  dynamics;  In  most  cases  the  analytical  work  has  not 
yet  provided  sequences  of  amino  acid  residues  and  Is  far  from  yielding  the  detailed 
picture  of  protein  structure  necessary  for  a  treatment  of  the  molecular  dynamics. 

For  this  reason  we  have  focussed  our  research  on  a  structural  feature  common  to  all 


2 


folded  or  helical  proteins#  the  amlde-I  hydrogen  bond  (-00-*  H-N-C=) .  The  dynamics 

of  this  system  are  complex,  due  to  electron  delocalization  on  the  -C-N-C=0  portion 

of  the  protein  backbone  and  the  nonlinear  effective  potential  at  the  proton.  Both 

phonon  spectra  and  lines  In  the  amide  sequence  of  acetanilide  exhibit  a  complex 

pattern  of  spectral  shifts  at  low  temperatures  (Johnston,  (7)).  The  pattern  of 
these  shifts  Indicates  the  presence  of  coupling  between  proton  vibrations,  acoustic 
phonons  and  excitations  of  the  C=0  and  C-N  bonds.  At  the  same  time  the  evidence 
Indicates  that  the  effective  potential  energy  of  the  proton  Is  a  double  well.  The 
classical  or  semlclasslcal  dynamics  of  such  nonlinear  systems  are  known  to  be  very 
sensitive  to  driving  fields  at  particular  frequencies.  Because  hydrogen  bonding  In 
crystalline  acetanilide  Is  a  good  qualitative  model  for  hydrogen  bonding  In 
proteins,  we  should  expect  to  see  spectral  features  In  substances  composed  of  alpha- 
helices  similar  to  those  found  In  acetanilide  crystals.  This  observation  would  have 
definite  Implications  for  the  dynamics  of  protein  molecules.  In  particular,  many  of 
the  normal  modes  of  a  protein  molecule  are  determined  by  parameters  arising  from  the 
hydrogen  bonding.  If  the  (-C=0** H-N-C=)  system  shows  high  sensitivity  to  driving 
electromagnetic  fields  In  the  millimeter  wave  range,  this  should  be  reflected  In  the 
dynamics  of  the  whole  molecule.  As  an  Initial  step,  we  have  chosen  to  Investigate  a 
simplified  model  proposed  by  Davydov  (8)  of  Interactions  between  amlde-I  excitations 
and  acoustic  phonons  In  the  protein  alpha-helix. 

A  major  thesis  Investigated  Is  that  the  principal  response  of  biological 
macromolecules  to  millimeter  wave  radiation  Involves  the  excitation  of  acoustic 
phonons  or  low  frequency  vibrations  with  characteristic  frequencies  less  than  6000 
GHz  (200  cm-1.)  The  biological  Importance  of  defining  the  properties  of  these  modes 
and  their  Interaction  with  electromagnetic  radiation  Is  stressed  by  Indications  of 
their  role  In  the  activity  of  enzymatic  proteins  (Petlcolas,  (9),  Northrup,  et  al  , 
(10))  and  In  directing  the  paths  of  hydrogen  exchange  In  proteins  (Sheridan,  et  al 
(ID).  As  a  result,  the  understanding  of  the  low  frequency  dynamics  of  proteins  and 
polypeptides  Is  the  subject  of  much  current  experimental  (Petlcolas,  (9))  and 
theoretical  research  (Northrup,  et  al.  (10)  Sheridan,  et  al,  (11),  McCammon  and 
Karplus.  (12)  Pain,  (13),  Go,  et  al,  (14),  Brooks  and  Karplus,  (15),  Davydov, 

(16, a, b),  Hyman,  et  al,  (17),  Scott,  (18),  Lomdahl,  (19)).  The  majority  of  the 
theoretical  Investigations  have  concentrated  on  analyzing  the  low  frequency  dynamics 
In  terms  of  such  weak  Interactions  as  hydrogen  bonding,  torsional  twisting,  and  van 
der  Waal's  forces.  Only  a  few  of  these  studies  Included  a  nonlinear  coupling 
between  the  low  energy  acoustic  phonons  and  the  amlde-I  vibrational  excltons. 

Part  II  of  this  report  covers  our  Investigation  of  the  dynamics  of  the  Davydov 
model  In  both  short  and  long  chain  alpha  helices  while  Part  III  summarizes  results 
of  molecular  orbital  calculations  of  the  excl ton- phonon  coupling  In  the  formamlde 
dimer.  Our  Investigation  of  these  two  aspects  of  the  Davydov  model  has  provided  a 
qualitative  understanding  of  oscillatory  modes  of  proteins  which  correspond  to  the 
millimeter  and  submil  1 imeter  regions  of  the  electromagnetic  spectrum. 

II.  INVESTIGATION  OF  THE  DYNAMICS  OF  THE  DAVYDOV  MODEL 

A.  BACKGROUND 

One  of  the  central  issues  of  bioenergetics  Is  how  energy  arising  from  chemical 
transformations  is  transferred  within  enzymes  and  other  protein  molecules  (Green, 
(20,21)).  The  amlde-I  vibration  (C*0  stretch)  has  been  proposed  as  a  possible 
vehicle  for  energy  propagation  In  biological  systems  (Green,  (21)).  A  linearized 
model  of  amlde-I  excitation  predicts  decay  to  thermal  vibrations  within  a  few 


picoseconds.  This  lifetime  seems  to  eliminate  propagation  of  amlde-I  excitation  as 
an  Intermediate  step  In  energy  transfer.  Most  chemical  processes  In  enzyme  systems 
Involve  conformational  changes  which  proceed  on  relatively  leisurely  timescales- 
microseconds  or  longer. 

Oavydov  (8)  has  proposed  a  potential  resolution  of  this  disparity.  He 
suggests  that  the  nonlinear  character  of  the  Interatomic  forces  can  stabilize  the 
amlde-I  excitations.  According  to  this  model#  propagation  of  an  amlde-I  excitation 
Is  analogous  to  an  exclton  In  a  one  dimensional  solid.  Excltons  are  coupled  to 
longitudinal  sound  waves  In  this  system.  This  coupling  of  phonon  and  exclton  fields 
produces  a  dynamically  stable  excitation  which  propagates  as  a  solitary  wave.  In 
essence#  the  phonons  provide  a  potential  well  that  prevents  exclton  dispersion. 

This  balance  between  dispersion  and  wave  steepening  Is  familiar  to  students  of 
solltonlc  phenomena.  In  the  continuum  approximation#  the  dynamics  of  the  solitary 
wave  are  described  by  the  nonlinear  Schroedlnger  equation#  which  has  solltonlc 
solutions.  Davydov's  analysis#  based  on  the  continuum  approximation#  demonstrates 
that  this  mechanism  could  stabilize  molecular  excltons  over  periods  of  milliseconds 
or  longer  (Davydov#  16a). 

Davydov#  Eremko  and  Serglenko  (22)  derived  a  phonon- exclton  model  which  more 
accurately  represents  the  three  dimensional  arrangement  of  hydrogen-bonded  chains 
along  the  protein.  Extensive  numerical  calculations  based  on  this  model  confirm 
that  Davydov's  continuum  approximation  Is  qualitatively  correct  for  alpha  helical 
chains  containing  several  hundred  peptide  residues.  All  models  exhibit  solitary 
waves#  which  form  at  some  threshold  value  for  the  coupling  between  excltons  and 
phonons.  Further#  the  solitary  waves  are  very  stable  and  closely  resemble  solltons. 

There  are  several  Issues  which  need  to  be  resolved  before  Davydov's  model  can 
be  rigorously  applied  to  proteins  In  a  living  cell.  These  Issues  fall  Into  four 
categories.  The  first  Issue  Is  the  range  of  chemical  or  dynamical  phenomena  to 
which  the  Davydov  model  can  be  applied  legitimately.  The  second  Issue  Is  whether  a 
model  related  to  the  Davydov  model  may  be  useful  for  structurally  complex  proteins. 
The  third  Issue  Is  the  relevance  of  the  parameters  used  In  previous  modeling 
efforts.  The  strength  of  the  exclton- phonon  coupling  Is  particularly  significant. 
The  fourth  Issue  Is  the  behavior  of  the  model  at  temperatures  commonly  found  in 
living  organisms. 

Davydov  originally  proposed  his  model  to  account  for  muscle  contraction 
(Davydov#  (8)).  Because  muscle  contraction  Is  thought  to  Involve  a  conformational 
change  In  a  long-chain  protein#  we  must  determine  whether  the  Davydov  model  allows 
for  a  substantial  portion  of  the  excltonlc  energy  to  be  transferred  to  phonons  or 
large-scale  nuclear  motions.  Unfortunately#  numerical  calculations  show  that  only  a 
small  portion  of  the  total  energy  resides  In  the  phonons  associated  with  the 
solitary  wave  (Scott#  (18)).  Even  though  the  Davydov  model  closes  the  gap  in 
timescales  between  dispersion  of  amlde-I  excitation  and  conformational  changes  in 
proteins#  the  gap  In  energy  remains  quite  wide.  Because  energy  transfer  from  an 
exclton  mode  to  a  highly  excited  phonon  mode  Involves  many  quanta#  the  probability 
of  energy  transfer  from  an  exclton  mode  to  a  highly  excited  phonon  mode  remains 
quite  small  (Ovchinnikov  and  Erlkhman#  (23)).  Nevertheless#  the  Davydov  model 
remains  a  reasonable  candidate  for  the  transfer  of  bond  energy. 

The  necessity  of  application  of  the  Davydov  model  to  a  wide  range  of  protein 
structures  arises  because  only  a  small  fraction  of  the  proteins  In  living  systems 
exist  as  long  alpha-helical  chains.  More  general  models  of  excltonlc  propagation 


are  necessary  for  the  highly  folded  forms  which  are  found  In  enzymes, 

mlcrof 11  aments,  and  chemical  receptors.  Lomdahl  (19)  has  proposed  one  such  model. 

This  model  reduces  to  the  Davydov  model  for  short  alpha  helical  segments. 


The  dynamics  of  the  Davydov  model  are  very  sensitive  to  the  effective  force 
represented  by  the  phonon-exclton  coupling.  Previous  numerical  studies  have  taken 
the  value  for  the  formamlde  dimer  calculated  by  Kuprlevlch  and  Kudrltskaya  (24)  as 
representative.  This  value  was  obtained  from  an  ab-initio  self-consistent  field 
computation.  The  authors  did  not  report  on  the  details,  which  Is  unfortunate 
because  such  calculations  are  quite  sensitive  to  choice  of  the  equilibrium  geometry, 
basis  sets,  or  configuration  Interaction  level.  Our  calculations  Indicate  the  value 
reported  by  Kuprlevlch  and  Kudrltskaya  may  underestimate  the  true  value  by  as  much 
as  a  factor  of  two.  We  conclude  that  an  accurate  assignment  of  the  coupling 
constant  remains  to  be  established  (See  part  III,  below). 

Another  related  problem  arises  In  the  extrapolation  of  experimental  results 
obtained  with  hydrogen-bonded  molecular  crystals,  such  as  acetanilide  (Carerl,  si 
il  ,  25))  to  support  the  Davydov  model  for  proteins.  The  Davydov  model  Includes 
only  excl ton- phonon  coupling  terms  which  are  linear  In  the  phonons.  Many  of  the 
models  of  propagation  of  excltons  In  molecular  crystals  Include  coupling  terms  which 
are  quadratic.  These  terms  are  necessary  to  account  for  phonon  frequency  changes 
which  may  be  present  In  addition  to  changes  In  the  equilibrium  positions  of  the 
molecules  (Munn  and  Sllbey,  (26)).  Spectrographlc  results  obtained  at  cryogenic 
temperatures  (Johnston,  (7,a,b)),  Indicate  that  the  phonon  frequency  in  acetanilide 
Is  shifted  by  a  change  In  the  position  of  the  proton  In  the  hydrogen  bonds  of 
crystalline  acetanilide. 

The  thermal 1z at Ion  problem  has  been  addressed  by  Davydov  (Davydov,  (27)).  By 
using  a  density  matrix  approach,  he  has  been  able  to  show  that  the  effective  phonon- 
exclton  coupling  constant  Is  decreased  by  the  presence  of  thermal  phonons  (Davydov, 
27,  28).  The  magnitude  of  this  effect  remains  to  be  computed.  In  fact,  there  Is  no 
conclusive  evidence,  either  theoretical  or  experimental,  which  would  Indicate 
whether  Davydov  solltons  exist  at  physiological  temperatures.  Our  results  Indicate 
that  thermal  effects  may  prevent  the  formation  of  Davydov  solltons  In  vivo  . 

B.  THEORETICAL  DERIVATION  OF  THE  MODEL 


Davydov’s  original  model  (Davydov,  (8))  Is  a  chain  of  N  molecules  along  the  z- 
axls.  The  equilibrium  positions  are  given  by: 


nR_  +  Un 
o  n 


(1) 


where  RQ  Is  the  equilibrium  distance  between  molecules  and  un  Is  the  displacement  of 
the  nth  molecule  from  equilibrium.  The  Hamiltonian  operator  for  excltons  (amlde-I 
bond  excitations  In  the  alpha  helix)  may  be  written: 

In  this  formula  £  Is  the  molecular  excitation  energy  and  J  is  the  energy  of 
resonant  Interaction  of  neighboring  molecules.  These  quantities  may  be  described  In 
terms  of  a  two-level  model  of  the  amlde-I  bond.  (Fedyanln  and  Yakushevlch, 

(29, a, b)).  In  this  representation  the  states  of  the  bond  are  denoted  0  and  f.  The 
excitation  energy  of  an  Isolated  molecule  Is  written  as  £  . 
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Then : 


t  -  ie+20l  D  -  <fO|Vn,Ml|fO>  -  <001  Vn,n*1|00> 
1  * 

We  may  also  write: 


where  d  Is  the  dipole  moment  of  the  molecules.  The  operators  Pp  and  P  +  are  the 
annhllatlon  and  creation  operators  satisfying  the  Pauli  commutation  relations: 
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The  phonon  Hamiltonian  has  the  form 
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where  pn  Is  the  momentum  conjugate  to  the  displacement  un: 

‘V  P.  J  (7> 

If  the  mass  of  the  molecule  Is  M  then 

w  *  Mv*2  =  (Mvk2)/(Rd2  (8) 

where  w  Is  the  elastic  coupling  constant  of  the  chain  and  v*1s  the  sound  velocity. 
The  dispersion  relation  for  phonons  Is  given  by 


flq  =2vJs1n(qR0/2)! 


(9) 


where  q  Is  the  wave  vector  and  Is  the  phonon  frequency.  The  Interaction  of 
excltons  with  phonons  Is  given  by 


Hint 
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where  "X  Is  the  exclton-phonon  coupling  parameter, 
model  Is 


(10) 

The  total  Hamiltonian  for  the 
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A  particular  solution  for  the  time-dependent  Schroedlnger  equation  associated  with 
Hso1  1s  91ven  by  the  following  wave  function: 


|fs(t)>  «  ItCn(t)  «xp{o(t)}P*|0> 
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(12) 
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In  this  expression  Cn(t)  Is  the  probability  amplitude  of  an  exclton  occupying  r-:*« 
n  at  time  t.  Conservation  of  energy  Implies  the  normalization  condition 


iicn(t)]2  *  1  . 

n 

The  operator 


(13) 


o(t)  =  i/f,I[£n(t)Pn  -  rn(t)unl  (14) 

Is  derived  from  a  coherent  state  representation  for  the  phonons.  This  follows  from 
the  relations: 


un  =  J(Tl/2MN&q)H  (bq  4  b-q ^  e*P 

+  (15) 

p  =  -1r(hf;  M/2N) 2  (b  -  b  )  exp  (inq; 
n  q  q  -q 

q  _ 

Substitution  of  these  Identities  Into  the  expression  for  <T(t)  gives  an  expression 

of  the  form 


eXP[1/t(Vq(t)bq  +  Z*q{tK]) 

where 

2c  ■  ^r,  exp  * 

q  r, 


(16) 


(17) 


i(ti'  M/2N)^2ir  exp  (inq) 

^  r 

The  properties  of  the  coherent  state  representation  Imply  that  the  functions 
TTn(t)  and  ^-_(t)  represent  the  average  momentum  and  displacement  respectively 
of  the  nth  molecule. 


Explicit  expressions  for  the  functions  /^n(t) *  7Tn(t >  and  Cn(t)  can  be  derived 
via  Hamilton's  canonical  equations  applied  to  the  energy  functional.  The  energy 
functional  Is  defined  as  follows: 


F  *  <Hsol>  *  <»#(t)|M!»s(t)>  .  (18) 

Me  may  calculate  F  from  the  definitions: 

F  *  ar|Cn(t)|2  -  JiiCJ+1(t)C  (t)  ♦  Cn+1(t)C*(t)] 
n  n 

+  (‘#<)i^(t)  +  w: [b  (t)  -  eri.1(t)]? 

h  In 

♦xjlcn(t)|2ient,(t)  -  6n.,(t)]  .  (19) 


The  canonical  Hamilton's  equations 


(20) 
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yield  the  identify 


Rlfn*l(t)  +  £n-l(t)  '  %•'>>  * 
(ft)!!C„+,(t)|2  -  |C„.,(t)|2, 


If  we  write  the  phonon  energy 
1 


w  =  4- 


2^n(t>  -  fe|e„(t)  -  tn.,(t)!2  , 


then  Schroedlnger 's  equation  Implies 

ttiCn(t)  3£. 

Tt  5C*(  t ) 

u  -  «>,un+1(t)  -en_,(t)i>  c„(t)  - 


(22) 


(23) 


Equations  (21)  and  (23)  describe  the  dynamics  of  a  single  chain. 


The  equations  for  a  single  chain  may  be  generalized  to  the  case  of  three 
Interacting  chains  as  In  a  protein  alpha  helix.  In  this  model,  resonance 
Interactions  may  occur  along  the  protein  helix.  Adding  these  Interactions  to  the 
model  yields  a  system  of  coupled  equations  In  which  the  dynamic  variables  a~e 
Indexed  by  spine  number  as  well  as  position  along  each  hydrogen-bonded  spine. 


Our  calculations  are  based  on  a  coupled  system  of  ordinary  differential 
equations,  obtained  as  outlined  In  the  previous  section.  Except  for  choice  of 
boundary  conditions,  our  procedures  are  formally  equivalent  to  those  of  Hyman, 
McLaughlin  and  Scott  (17).  The  ends  of  the  helix  are  not  mechanically  constrained. 
The  complete  model  Is  given  by  the  following  system  of  equations 


idC 


n,i  = 


r 


n+1  ,i 


-  0-'>56(Vl,itVl,1>  +  0-092lCn,W,Cn,i.l> 


♦  °-372  Oo10x)Ki-«")  -  sM>cn,i 

N 

*  <>  -0  11 
^n+l,i  "  6n,i^  +  ^  4o^Cn-l ,i ^£n,i  ’ 


8 


(24) 


d2a 

dt,2  U  y  1  *n+l  ,i  “n ,i 
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♦  <K>  <V,1  -  Vl.il  ♦  0.132(10  >)  x 


<’•*;>  «’*«:»  >1^1, i  i2-  Kno.ii2  > 


«>!,/  *  !Co.i:‘l  '  ^iCn-l,ii‘  *  lCn,i 


+  C*  .(C  ,  .  -  C  ,  .)  +  (C*  ,  .  -  C*  ,  .)c  .} 

n+l,v  n+l,i  n-l,v  '  r+1  ,i  n-1  ,v  n,i 


We  have  normalized  the  time  variable  In  these  equations  by  setting  t1  B  (w/M)*^  t. 
The  quantities  S0n,  £Nn  account  for  the  boundary  conditions  at  the  n  *  0  and  n  =  K 
ends  of  the  helix  respectively. 

For  the  computations  on  short  alpha  helices,  one  quantum  of  amlde-I  excitation 
energy  was  assumed  to  be  located  at  the  (0,1)  position  at  the  Initiation  of  the 
simulation.  No  provision  was  made  for  random  motions  due  to  thermal  disorder  during 
the  course  of  the  computations.  As  In  the  work  of  Scott  (18),  the  only  free 
parameters  of  the  model  are  the  number  of  peptide  residues  In  the  alpha  helix  and 
the  nonlinear  coupling  coefficient  between  excltpps  and  phonpps.  During  the  course 
of  the  computations,  we  used  values  of  0.2  X  10”  *  0.4  X  10“  »  and  0.6  X  10”*" 
Newtons  for  chain  lengths  of  18,  24,  27,  30,  and  60  peptide  residues. 


For  the  calculations  on  the  long  alpha  helices  one,  two  and  three  quanta 
Initializations  for  the  exclton  energy  were  employed.  Chain  length  was  200 
residues,  coupling  constant  ranged  from  0.5  X  10”10  Newtons  to  3.0  x  10”*°  Newtons. 
Calculations  were  carried  out  for  the  chain  Initially  at  rest  and  for  phonon  motion 
simulating  a  temperature  of  300°  Kelvin. 


The  resulting  systems  of  differential  equations  were  solved  on  a  VAX  11/780 
using  double  precision  arithmetic.  All  Integrations  were  performed  using  a  third- 
order  Adams-Bashford-Moulton  predictor-corrector  method  (Shamplne  and  Gordon,  (30)). 
The  code  was  set  to  maintain  a  relative  error  of  0.000001  at  each  time  step  during 
the  computation.  Computation  of  the  total  probability  of  bond  occupancy 


served  as  an  independent  check  of  accuracy.  In  every  computation  P  -  1  was 
maintained  at  a  value  less  than  0.00001. 
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C.  RESULTS 


1.  SHORT  ALPHA  HELICES 

The  total  phonon  energy  was  plotted  as  a  function  of  time  in  order  to  examine 
the  energy  exchange  between  excltons  and  phonons.  Examination  of  the  plots  for  the 
fifteen  combinations  of  chain  length  and  coupling  constants  described  above  shows 
that  several  time  scales  are  present: 

a)  Short  period  oscillations  (about  3-6  terahertz).  These  may  be  seen  as 
transient  bursts  which  recur  at  periodic  Intervals. 

b)  The  time  Interval  between  the  appearance  of  the  bursts  Is  roughly 
proportional  to  the  chain  length.  The  duration  for  each  burst  appears  to  be 
independent  of  chain  length.  Figures  1-3  show  plots  of  total  phonon  energy  for 
chain  lengths  of  18*  27  and  36  residues.  The  high  frequency  oscillations  recur 
at  3  X  10  12  seconds  for  chain  length  18*  4.5  X  10“1Z  seconds  for  chain  length 
27  and  6  X  10"1Z  seconds  for  chain  length  36. 

c)  Intermediate  period  fluctuations  (500  gigahertz  to  1  terahertz).  The  time 
period  for  these  fluctuations  is  proportional  to  chain  length.  Intermediate 
period  fluctuations  are  marked  in  Figures  2  and  3. 

d)  Long  period  fluctuations  (100  gigahertz).  Long  period  fluctuations  are 
seen  for  chain  lengths  of  27  and  36  residues.  The  clearest  example  may  be  seer. 
In  Figure  2  for  a  chain  length  of  27  residues,  with  coupling  coefficient  C.4  X 
10“10  Newtons.  This  oscillation  persists  for  many  cycles  In  the  original 
record  from  which  the  figure  Is  extracted. 

Reordered  according  to  frequency,  the  four  separate  timescales  observed  are: 

a)  Fast  oscillations,  3  to  6  terahertz,  occurlng  In  bursts. 

b)  Intermediate  fluctuations,  500  to  1000  gigahertz. 

c)  Recurrence  of  bursts,  200  to  500  gigahertz. 

d)  Long  period  fluctuations  below  200  gigahertz. 

Examination  of  the  modules  of  the  Fourier  transform  of  the  total  phonon  energy 
confirms  these  observations.  As  the  chain  length  Is  Increased,  the  short  period 
oscillations  become  less  significant  and  the  long  period  oscillations  more 
significant.  This  trend  Is  shown  in  Figures  4  and  5.  The  division  between  the 
timescales  Is  seen  most  clearly  for  a  chain  length  of  27  residues. 

Phonon  Spectra 

Because  phonon  energy  In  an  actual  protein  alpha  helix  Is  not  easily  subject  to 
direct  experimental  observation,  we  also  computed  the  7T(t)  momentum  expectation 
function  for  selected  residues  along  the  chain.  The  modulus  of  the  Fourier 
transform  of TT(t)  gives  the  phonon  spectrum,  as  may  be  seen  from  equation 
(15)  and  the  relation TTn(t)  *  <Ys(t)|pn|  Y*lt)>  (See  Haaken,  31).  In  order  to 
analyze  the  collective  modes,  we  evaluated  the  sum  of  momenta  for  residues  1,  2  and 
3  and  the  sum  of  momenta  for  residues  [n/2],  [n/23  -  1  and  [n/23  -  1  and  Cn/23  -  2. 
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Figure  1. 


Phonon 


Coupling  coefficients: 

0.6  X  10-10  Newtons. 

Bursts  of  6  Thz  oscillations  may  be  seen  at  time  Intervals  of 
approximately  3.2  X  10“12  seconds.  Bursts  are  Indicated  by  arrows. 


n  energy  for  an  18  residue  chain. 
.2  X  10"10,  b)  0.4  X  10’10,  c) 
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Figure  2.  Phonon  energy  for  a  27  residue  chain.  Excl ton- phonon  coupling 

coefficients  are  as  in  Figure  1).  Bursts  of  6  Thz  oscillations  may 
be  seen  at  time  Intervals  of  approximately  4.8  X  10”1Z  seconds. 
Location  of  bursts  Indicated  by  arrows  above  traces.  Periodicity  of 
Intermediate-time  scale  fluctuations  Is  Indicated  by  arrows  below 
the  traces.  A  long-period  fluctuation  Is  Indicated  on  b). 
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Figure  4.  Fourier  transforms  of  total  phonon  energy.  Phonon  energy  for  chains 
of  a)  18  residues*  b)  27  residues*  and  c)  36  residues.  Nonlinear 
coupling  coefficient  Is  0.2  X  10”10  Newtons.  Note  the  trend  toward 
lower  frequencies  for  longer  chains. 
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Figure  5.  Fourier  transforms  of  total  phonon  energy  data.  Phonon  energy 
for  chains  of  a)  18.  b)  27  and  cl  36  residues.  Nonlinear 
coupling  coefficient  Is  0.6  X  10“iO  Newtons.  Frequency  scale  is 
as  In  Figure  4. 
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Samples  of  the  momenta  plots  are  shown  In  Figure  6.  We  also  computed  the 
modulus  of  the  Fourier  transform  of  each  momentum  plot.  Plots  of  amplitude  versus 
frequency  are  shown  In  Figures  7,  8  and  9.  Each  of  the  plots  consists  of  a  series 
of  well-defined  peaks.  These  peaks  correspond  to  the  normal  modes  of  the 
corresponding  linear  chains*  and  follow  the  same  dispersion  relation  as  given  In 
equation  (9). 

2.  LONG  ALPHA  HELICES 


Our  finding  on  long  alpha  helices  are  summarized  In  Figures  10-17.  (Each  of 
these  plots  required  more  than  10  hours  of  CF’U  time  to  produce.)  As  may  be  seen 
from  Figures  10*  11  and  12*  representing  evolution  of  probability  of  exclton 
occupation  subsequent  to  Initial  excitations  of  one.  two  and  three  quanta  of  amide- I 
excitation  energy  placed  In  the  end  peptide  groups  of  chains  at  0°Kelv1n,  sollton 
formation  Is  a  function  of  excitation  level.  Even  In  the  case  of  two  or  three 
quanta,  there  Is  a  significant  probability  that  a  fast-moving  dispersive  wave  will 
form.  This  may  be  seen  by  the  sharply-defined  ripples  which  form  upon  reflection  of 
the  wave  from  the  end  of  the  chain.  Another  branch  of  the  dispersion  of 
excitations*  which  Is  Indicated  by  waves  similar  to  but  moving  more  slowly  than  the 
main  solitary  wave.  Note  that  the  coupling  constant  In  this  case  Is  0.5  X  ltT10 
Newtons*  a  value  which  Is  reported  by  Scott  (18)  to  be  In  the  middle  of  the  range  of 
coupling  constants  at  which  solitary  waves  form. 


The  second  series  of  Figures,  13*  14  and  15.  Illustrate  the  effects  of  high 
coupl Ing  constants.  Figure  13,  which  represents  evolution  of  exclton  probability  In 
a  chain  at  0  degrees  Kelvin  with  a  coupling  coefficient  of  1.5  x  10“10Nt,  shows  that 
a  slowly-moving  solitary  wave  may  form  at  a  lower  excitation  level  In  the  case  of 
stronger  coupling.  Figure  14,  also  representing  exclton  evolution  at  0  degrees 
Kelvin  but  at  a  coupling  constant  of  3.0  x  10“1DNt,  shows  that  at  very  high  levels 
of  excitation,  a  randomly-moving  exclton  forms.  This  Is  familiar  to  solid-state 
physicists.  In  that  phonons  scattering  from  the  exclton  push  It  from  site  to  site. 
Figure  15.  representing  evolution  subsequent  to  a  two  quanta  excitation,  shows 
eventual  trapping  of  the  exclton  at  the  end  of  the  chain.  A  similar  phenomenon  has 
been  reported  by  Scott  (18),  with  the  major  difference  that  the  exclton  did  not 
Initially  leave  the  end  of  the  chain. 


Figures  16  and  17  illustrate  the  effects  of  thermal Izatl on.  Figure  16,  shows 
the  evolution  subsequent  to  a  two  quanta  excitation  In  a  chain  with  thermal  phonon 
motion  representing  a  temperature  of  238°  Kelvin.  The  excltons  disperse  rapidly, 
within  20  picoseconds.  Figure  17  shows  exclton  evolution  subsequent  to  three  quanta 
excitation  In  a  chain  with  a  phonon  temperature  of  300  degrees.  Dispersion  Is  also 
rapid  In  this  case. 


D.  DISCUSSION 
1.  SHORT  ALPHA  HELICES 

We  believe  that  the  full  range  of  phenomena  observed  during  the  course  of  our 
modeling  studies  requires  a  more  detailed  analysis  than  the  passage  to  the  continuum 
limit  which  Is  common  to  the  literature  on  the  Davydov  model  (Davydov,  8,  16,  27, 

28,  Scott,  18).  A  more  accurate  treatment  should  follow  the  spirit  of  Toda's 
analysis  of  nonlinear  lattices.  (Toda,  32). 
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Figure  6.  Velocity  of  displacement  for  groups  at  mid-section  of  alpha  helices 
of  a)  18'  b)  30,  c)  36  and  d)  60  residues.  Nonlinear  coupling 
coefficient  Is  0.2  X  10“10  Newtons. 
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Figure  7.  Phonon  spectra.  Chain  length*  24  residues.  Site  Is  mldchaln. 

Nonlinear  coupling  coefficients  are  a)  0.2*  b>  0.4*  c)  0.6 
XIO"10  Newtons.  Note  that  relative  amplitudes  and  positions  of 
peaks  remain  the  same.  Horizontal  scale  In  terahertz.  Vertical 
scale  Is  arbitrary. 
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Figure  8.  Phonon  spectra.  Same  scales  as  In  Figure  7).  Chain  lengths  and 
sites  are  a)  18  residues*  end*  b)  18  residues*  middle*  c)  36* 
middle*  d)  36*  end.  Coupling  coefficient  Is  0.4  x  10“iONt.  No 
that  peaks  are  In  the  correct  harmonic  relationship. 
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Figure  9.  Phonon  spectra.  Same  scales  as  In  Figure  7).  Chain  lengths  and 

sites  are  a)  30*  end*  b)  30*  middle*  c)  60*  end*  d)  60*  middle 
Coupling  coefficient  Is  0.4  x  10“10Nt. 
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Figure  10.  Plot  of  probability  of  exclton  occupation.  Horizontal  axis  Is 
position  along  a  chain  of  200  residues.  Each  trace  represents  a 
time  Increment  of  approximately  0.4  picoseconds.  Phonon- exclton 
coupling  parameter  Is  0.5  x  10  10Nt.  Note  that  most  of  the 
excl tonic  energy  propagates  as  a  dispersive  wave. 
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Figure  11.  Plot  of  probability  of  exclton  occupation.  Coupling  constant  of 
0.5  X  10_1®Newtons*  two  quanta  Initialization.  Note  the  presence 
of  both  a  solitary  wave  and  a  faster-moving  dispersive  wave. 
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,  Figure  12.  Evolution  of  exclton  occupation  probability.  Coupling  constant  of 

|  0.5  X  10”10  Newtons#  three  quanta  Initialization. 
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Figure  14.  Evolution  of  exclton  occupation  probability.  Coupling  constant  of 
3.0  X  10"10  Newtons.  One  quantum  Initialization.  Note  that  a 
solitary  wave  forms  In  this  case. 


25 


RESIDUE  NUMBER 


Evolution  of  exclton  occupation  probability.  Coupling  constant  of 
3.0  X  10"10  Newtons.  Two  quanta  Initialization.  Note  that  a  self 
trapped  exclton  Is  attracted  to  chain  boundary. 
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Figure  16.  Evolution  of  exclton  occupation  probability.  Coupling  constant  of 

C  Y  UoulAfic  Th.AA  Aii.nl.  Inltlil  l.illnn  Phftnnn 


0.5  X  10"iU  Newtons.  Three  quanta  Initialization.  Phonon 


temperature  of  300  Kelvin.  Note  rapid  dispersion  of  excltons. 
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Figure  17. 
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Evolution  of  exclton  occupation  probability.  Coupling  constant  of 
0.5  x  10”10  Newtons.  Two  quanta  Initialization.  Phonon  temperature 
of  238°  Kelvin.  Excltons  disperse  rapidly. 
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Equations  (24)  may  be  simplified  by  neglecting  coupling  terms  of  small 
magnitude.  In  the  first  equation,  this  includes  the  observation  that  only  the 
hydrogen  bond  on  the  OO  side  of  the  peptide  unit  will  influence  the  amlde-I  energy: 


(25) 


idC. 
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dt’ 


d2£ 


dt 


.c 


■  C.372O0,0x)l(r„4lii-  6„,,)Cni,  *  U„,r 
-  0.0S8(Cn+1>1  *  *  0-092(C„iftl  *C„i10) 
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For  the  sake  of  simplicity,  we  have  not  Included  the  end  conditions  In  these 
equations.  These  equations,  which  have  several  features  In  common  with  Flaschka's 
equations  for  the  Toda  lattice  (Flaschka,  33),  may  be  treated  as  matrix  equations. 
To  begin  the  analysis,  note  that  to  the  first  order  the  equations  for  the  excltons 
and  phonons  are  uncoupled.  The  equations  for  the  excltons  may  be  written  as  a 
system  of  linear  autonomous  equations: 

1/dA/dt'  *  KA,  (26) 

with  solution  A  *  exp(1K)A0. 


The  phonon  displacement  Is  given  In  turn  by  an  approximate  equation: 

1d2£/dt'2  *  K&+  F,  (27) 

where  (5  n> j  Is  given  by  1Bn  ,+  B_j)  and  F  Is  the  nonlinear  term  given  by  the 
excltons/  The  solution  of  (27)  Is  given  by  the  formula 

(B  (t)  *  Jo^xpC-mit-TJDFITJd/:  (28) 

This  scheme  may  be  continued  self-conslstently  to  Improve  the  estimate  of  the 
evolution  of  the  exclton  probability  amplitude,  given  by  equation  (27),  but  equation 
(28)  accounts  for  the  major  features  we  have  noted  In  the  results  section.  For 
example,  the  average  lifetime  of  an  exclton  at  a  site  of  the  chain  Is  about  1.6 
picoseconds.  An  exclton  moving  along  the  chain  will  excite  1.6  phonons  most 
strongly  when  It  Is  near  one  of  the  ends  of  the  chain.  This  Is  the  origin  of  bursts 
In  the  phonon  mode  at  about  6  THz.  Other  fluctuations  are  due  to  resonances  between 
phonons  and  the  periodicity  of  the  excltons.  The  long  period  observed  for  the  chain 
length  of  27  residues  Is  due  to  transfer  between  even  and  odd-numbered  spines  of  the 
alpha  helix. 

Because  the  calculated  phonon  spectra  do  not  reflect  the  frequency  shifts 
observed  In  acetanilide  crystals,  we  suggest  that  the  Interaction  term  of  the 
Hamiltonian  should  Include  a  quadratic  term.  Writing 
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(30) 


would  yield  a  new  expression  for  the  energy  functional 


F‘  -  F  +  teic/I£n+l(t)  •  Vl(t)  ] 

n 

Equation  (23)  would  be  modified  by  addition  of  a  term 


f  A*/16 


n+2 
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£n;t>J 


Er  -  Vl(t)]} 


Similar  modifications  may  be  made  In  the  three-spine  model. 
2.  LONG  ALPHA  HELICES 
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Mathematical  analysis  of  the  complicated  dispersion  observed  for  a  single  quantum 
excitation  Is  not  yet  complete.  The  reason  for  sollton  formation  In  the  case  of  two 
and  three  quanta  excitation  may  be  seen  from  the  fact  that  phonon  formation  goes  as 
the  square  of  the  number  of  excltons  and  thus  the  "feed-back"  term  In  the  exclton 
evolution  equation  (25)  Increases  as  the  cube  of  the  number  of  excltons.  The 
qualitative  change  In  exclton  dynamics  observed  when  the  coupling  constant  reaches 
3.0  X10”*°Nt  is  also  due  to  a  qualitative  change  In  the  "feed-back"  term.  In  this 
case,  the  coefficient  of  the  feedback  term  exceeds  unity. 


Ill  CALCULATION  OF  THE  EXCITON-PHONON  COUPLING  OONSTANT 


A.  BACKGROUND 

The  calculation  of  the  non-linear  coupling  between  the  vibrational  excltons 
localized  In  amide- I  bonds  and  the  phonons  In  the  hydrogen-bonded  "spines"  of  a 
polypeptide  alpha-helix  Is  Important  to  the  understanding  of  the  propagation  of 
undispersed  amlde-I  vibrational  excltons  along  these  "spines".  An  undispersed 
exclton  which  propagates  through  an  Interaction  with  phonons  Is  also  defined  as  a 
Davydov  sollton.  In  a  numerical  analysis  of  a  polypeptide  alpha-helix  of  finite 
length,  Scott  showed  that  a  Davydov  sollton  can  be  sustained  when  the  exclton-phonor 
coupling  parameter  exceeds  a  minimum  magnitude  of  0.35  X  10“iONt(18).  The 
question  of  the  validity  of  this  threshold  Is  addressed  In  the  calculations 
described  below. 

According  to  Kuprlevlch  and  Kudrltskaya,  the  excl ton-phonon  Interaction 
parameter  consists  of  a  resonance  component  (X2)  and  two  non-resonance  components 

(xl+'x1‘)(24).  The  X2  term  can  be  obtained  from  the  known  magnitude  of  the 
resonance  Interaction  between  the  translstlon  dipoles  associated  with  the  amlde-I 
vibration  In  alpha-helical  polypeptides  (44).  In  this  manner.  Scott  derived  a  value 
of  10-12  Newtons  for  X_ (18).  The  non-resonance  components  are  best  determined  from 
a  theoretical  treatment  of  excl ton- phonon  coupling  at  the  microscopic  level.  This 
treatment  Involves  the  use  of  the  quantum  mechanical  ab  Initio  self-consistent  field 
(SCF)  molecular  orbital  (MO)  procedure.  The  excl ton-phonon  coupling  terms  and  the 
application  of  ab  Initio  SCF  MO  theory  to  the  calculation  of  these  terms  are 
discussed  In  greater  detail  In  Section  III  B.  Section  III  C  presents  our 
preliminary  calculation  of  Xj*  and  that  of  Kuprlevlch  and  Kudrltskaya  (24)  for  the 
simplest  hydrogen-bonded  polypeptide  -  the  formamlde  dimer.  Conclusions  are 
presented  In  Section  III  D. 

The  Hamiltonian  for  the  Interaction  of  vibrational  excltons  with  acoustical 
phonons  In  a  discrete  chain  of  N  molecular  units  forming  a  regular  one-dimensional 
lattice  was  derived  by  Kuprlevlch  and  Kudrltskaya  (24,  34),  In  a  form  differing 
slightly  than  that  reported  by  Davydov. 

•W  •  l2.ipSp.«i+<w  v’«f  Vi« 

w  .  .  (32) 

4  2l!=2<PmVl4fLlP.»<2'»»-Vl> 

In  this  expression,  P_+  and  Pm  are  the  creation  and  annihilation  operators  of 
vibrational  excitation  In  the  unit  m,um  Is  the  displacement  of  this  unit  from  Its 

equilibrium  position,  X.+  and  X^“  are  the  parameters  defining  the  non-resonance 
Interaction  of  a  unit  with  Its  neighbors,  and  X2  Is  the  parameter  defining  the 
resonance  1nteract1on_for  two  neighboring  units  In  a  chain.  Note  that  the  terms  In 
the  first  sum  with  Xj"  at  m=l  and  with  X  j  at  m»N  should  be  omitted,  because  both 

terminal  molecular  units  have  only  one  neighbor. 
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The  explicit  formulations  of  X,4,  *  end  X2  ^or  *  dimeric  group  consisting  of 

peptide  groups  m  and  m+1  are  as  follows  (24,  34): 

V  ■  (33) 

X,-  «  (V2.HI)‘E„Bl/)Vl!in,»(l  (34) 

X2  -  <*/2.)>(}1E,„,wl/j^Qntl>/iR|„>mtl  (35) 

where  w  Is  the  Intramolecular  amlde-I  vibrational  frequency,  FL  ,^1  Is  the  distance 
or  hydrogen-bond  length  between  groups  m  and  m+1  (See  Figure  lea),  Is  the 

Interaction  energy  (or  total  energy  In  the  Born-Oppenhelmer  approximation (35) )  of 
groups  m  and  m+1,  and  Q.^  Is  the  normal  mode  coordinate  for  the  amlde-I  vibrations  In 
groups  1*m  and  m+1.  Since  the  second  derivatives  of  the  total  energy  of  the  dimer 
with  respect  to  Its  vibrational  normal  coordinates  are  proportional  to  the  harmonic 
force  constants  for  these  modes,  the  coupling  parameters  In  Eqs.  (33)  -  (35)  can  be 
viewed  as  the  first  derivatives  of  the  harmonic  force  constants  for  the  amlde-I 
vibrations  In  the  dimer  with  respect  to  the  Intermolecular  distance  R_  Thus, 

the  acoustical  phonons  of  the  alpha-helix  which  are  defined  principally  by  changes 
In  R,,,^  d.e.,  changes  In  the  pitch  of  the  alpha-helix  (24,  34))  are  coupled  to 

the  amlae-I  vibrational  excltons. 

As  stated  In  the  Introduction,  a  value  of  10“12  Nt  for  X2  can  be  derived  from 
the  known  magnitude  of  the  resonance  Interaction  between  the  amlde-I  vibrational 
transition  dipoles  (18).  Kuprlevlch  and  Kudrltskaya  argue  that  Xj~  is  negligible 

with  respect  to  X,4  (24).  The  principal  reason  Is  that  the  dominant  contribution  to 
the  amlde-I  vibrational  normal  mode  In  a  peptide  group  Is  the  OO  bond  stretching 
motion  (24,  36,  37).  As  a  result,  the  amlde-I  vibration  of  group  m  will  be  much 
more  perturbed  than  that  of  group  nr+1  by  changes  In  This  observation  Is 

clearly  seen  In  Figure  18a.  The  major  objective  of  our  analysis  of  the  exclton- 
phonon  coupling  In  polypeptides  will  therefore  be  the  calculation  of  Xj4. 

B.  QUANTUM-MECHANICAL  CALCULATION  OF  THE  NON-RESONANCE  EXCITON-PHONON  INTERACTION 

PARAMETER  Xj. 

The  calculation  of  the  Interaction  energy  of  the  two  peptide  groups  (Em  ^j) 
as  a  function  of  the  Intra-  and  Inter-molecular  coordinates  Is  fundamental  to  the 
determination  of  Xj4.  In  the  Born-Oppenhelmer  or  adiabatic  approximation  (35),  the 
Interaction  energy  Is  equivalent  to  the  total  energy  (E)  for  the  system,  which  Is 
calculable  using  quantum  mechanical  ab  Initio  self-consistent  field  (SCF)  molecular 
orbital  (MO)  procedures.  The  details  of  these  procedures  and  their  application  to 
the  study  of  the  electronic  structure  of  atoms  and  molecules  are  well-covered  by 
Schaefer  In  his  book  (38)  and  the  references  therein. 

The  mapping  of  the  E  surface  for  a  dipeptide  system  as  a  function  of  Its  tntra- 
and  Inter-molecular  coordinates  Is  a  formidable  task.  However,  as  outlined  by 
Kuprlevlch  and  Kudrltskaya  (24),  one  can  obtain  a  crude  estimate  of  X.4  by  the 
following  procedure.  1)  Investigate  the  simplest  system  that  contains  two  peptide 
groups  connected  by  a  hydrogen  bond,  l.e.  the  formamlde  dimer  shown  1h  Figure  18b. 
2)  Determine  an  equilibrium  geometry  for  the  dimer.  3)  Approximate  the  amlde-I 
vibrational  normal  mode  as  only  a  C»0  bond  stretching  motion.  As  a  result,  the 
In  Eq  (33)  becomes  T7w  (24),  where  |Ms  the  change  In  C*0  bond  length  away  from 
Its  equilibrium  value  (positive/negative  J  signifies  an  expansion/contraction  of  the 
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Figure  18 

Figure  18 
Figure  18 


(IBi 


* 

/ 


* 


(a) .  Structure  of  a  chain  of  hydrogen-bonded  peptide  groups. 

Unspecified  molecular  fragments  are  designated  R  in  the 
figure.  Note  that  the  Integer  m  represents  one  peptide  group. 
The  hydrogen  bond  length  between  groups  m  and  m+1  is  Indicated 

*m,  «*!' 

(b) .  Structure  of  the  formamlde  dimer  and  numbering  of  the  atoms 

used  In  the  calculations. 

(c) .  Description  of  the  vibrational  normal  mode  In  the  formamlde 

monomer  which  contains  the  stretching  motion  of  the  OO  bond. 
This  figure  Is  taken  from  Ref  37,  and  the  direction  and  length 
of  the  arrows  represent  the  activity  of  the  given  atoms  In  the 
normal  mode. 
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C*0  bond),  K  Is  the  harmonic  force  constant  for  the  OO  bond  stretching  motion,  and 
»1s  the  amide- I  frequency.  Eq.  (33)  then  takes  the  form 

Xx+*  HW2K)  (36) 

•here  K  *  3*E/)5^.  4)  Calculate  E  at  three  different  values  of  5  and  estimate 
K  0*E/Jjr  by  d(dE)/(8f)  .  5)  Calculate  K  at  different  values  of  R^,  ^  and  + 

estimate  JK/i  IL  ^2  by  AK/&  1°  summary,  a  crude  determination  of  Xj 

for  a  polypeptide  chain  can  be  made  by  calculating  the  total  energy  of  the  formamlde 
dimer  In  different  geometries  and  then  taking  the  appropriate  numerical  derivatives. 

Both  Kuprlevlch  and  Kudrltskaya* s  (24)  and  our  ab  Initio  SCF  MO  calculations  of 
the  total  energy  of  the  formamlde  dimer  were  performed  with  the  minimal  atomic 
orbital  basis  set  composed  of  36  Slater-type  orbitals  (STO),  each  being  approximated 
by  3  Gaussian  functions  (3G)  (39).  It  Is  stressed  that  the  STO-3G  basis  set  Is  good 
only  for  qualitative  analyses  of  the  geometries  of  amide  compounds  (37,  38).  For 
example,  an  energy-minimized  geometry  of  the  formamlde  monomer  calculated  using  the 
STO-3G  basis  set  was  not  In  complete  agreement  with  the  experimentally  determined 
geometry  for  the  molecule  (37).  The  major  differences  were  that  the  calculated 
geometry  was  slightly  non-planar,  while  the  observed  one  was  planar,  and  that  the 
calculated  C-N  bond  length  was  0.1a  longer  than  the  observed  one  (see  Table  I  and 
Ref.  (37)).  Although  the  ST0-3G  basis  set  has  Its  limitations  with  regard  to  the 
accurate  study  of  molecular  properties,  SCF  convergence  was  achieved  for  all  6  of 
geometries  of  the  formamlde  dimer,  and  the  cost  per  calculation  was  about  12  minutes 
of  CPU  time  on  the  VAX  11/780.  Thus,  the  STO-3G  basis  set  allows  one  to  make  a 
qualitative  estimate  of  X2+  with  the  least  amount  of  computational  effort. 

The  specific  differences  between  Kuprlevlch  and  Kudrltskaya's  (24)  and  our 
approach  to  the  calculation  of  X2+  are  the  determination  of  an  equilibrium  geometry 
for  the  dimer  and  the  manner  In  which  J  was  changed.  The  equilibrium  geometry  of 
Kuprlevlch  and  Kudrltskaya  (24)  was  found  as  follows.  The  relative  arrangement  of 
nuclei  In  each  of  the  dimer  molecules  was  assumed  to  be  the  same  as  observed  In  a 
free  formamlde  molecule  (see  Table  I).  The  mutual  orientation  of  the  two  molecules 
was  constrained  so  that  the  atoms  C2*0i — H4-N2  on  *  straight  line  (see  Figure 
18b).  The  equilibrium  Intermolecul ar  distance  was  taken  as  that  corresponding  to 
the  minimum  total  energy  along  the  C2*0j — Hj-N2 

Our  equilibrium  geometry  was  derived  from  an  energy  minimization  of  the 
complete  geometry  of  the  dimer.  This  calculation  took  approximately  10  hours  of  CPU 
time  and  was  accomplished  utilizing  energy  gradient  and  optimization  routines.  The 
only  constraint  on  the  geometry  optimization  was  that  the  dimer  had  to  be  planar. 
This  constraint  was  Implemented  because  the  observed  geometry  of  the  formamlde 
monomer  Is  planar  (see  Table  I  and  Ref.  (40)). 

The  change  In  }  should  correspond  to  the  motion  of  the  carbon  and  oxygen  atoms 
In  the  C2*02  bond  stretching  motion.  Kuprlevlch  and  Kudrltskaya  (24)  chose  to  move 
only  the  0^  atom  In  Increments  of  -0.02  Bohr,  so  that  E  was  calculated  at  the  points 
5*0.0,  -0.02,  and  -0.04  Bohr.  The  reason  for^the  choice  of  two  negative  values  of  Jf 
Is  that  the  observed  C*0  bond  length  of  1.243A  (41)  used  as  their  soul librium  value 
Is  longer  than  the  one  calculated  using  the  ST0-3G  basis  set  (1.219A)  (see  Table  I). 
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We  elected  to  move  both  the  Oj  and  atoms  In  Increments  of  0.01  Bohr,  so  that 
E  was  calculated  at  JB0.0*  0.02*  and  0.04  Bohr.  Note  that  the  molecular  fragment 
attached  to  the  Ci  atom  (atoms  Nj,  Hj  *  H2*  and  Ho  In  Figure  18b)  was  also  moved 
rigidly  with  the  atom.  Our  motivation  for  moving  both  the  Oi  and  atoms  was 
that  this  motion  Is  In  better  agreement  with  the  amlde-I  vibrational  normal  mode  for 
the  formamlde  monomer  shown  In  Figure  18c. 

C.  RESULTS 

Calculated  and  observed  equilibrium  Internal  coordinates  and  ground  state 
dipole  moments  for  the  formamlde  dimer  and  monomer  are  given  In  Talkie  I.  The 
calculated  ground  electronic  state  properties  used  to  determine  for  the 
formamlde  dimer  are  presented  In  Table  II.  Table  III  lists  Kuprlevlch  and 
Kudrltskaya's  (24)  and  our  calculated  values  of  X,+  for  different  harmonic  force 
constants  for  the  C=0  bond  stretch  vibration  In  the  formamlde  dimer  and  monomer. 

Equilibrium  Geometry 

No  experimental  data  are  available  on  the  geometry  of  the  formamlde  dimer  shown 
In  Figure  18b.  In  order  to  estimate  the  accuracy  of  our  energy-mi nlmlzed  geometry 
of  the  formamlde  dimer  calculated  using  the  STO-3G  basis  set*  we  calculated  an 
energy-minimized  geometry  for  the  formamlde  monomer,  on  which  there  Is  experimental 
data  (40).  As  In  the  case  of  the  dimer,  the  monomer  was  constrained  to  be  planar 
(designated  ST0-3G  planar  In  Table  I).  We  see  In  Table  I  that*  except  for  the 
longer  Ci-Ni  bond  length,  our  calculated  geometry  does  not  differ  significantly  fror 
the  observed  geometry  given  In  the  rightmost  column  of  Table  I,  or  from  the  geometry 
calculated  using  the  ab  Initio  SCF  M0  procedure  with  the  more  elaborate  4-31G  split 
valence  basis  set  (42).  We  conclude  that  our  calculated  equilibrium  geometry  for 
the  formamlde  dimer  Is  reasonable  with  the  exception  of  the  Cj-N,  and  C2-N2  bond 
lengths.  Most  Importantly*  the  calculated  bond  length  of  1.218  A  In  the 

formamlde  monomer  Is  In  very  good  agreement  with  observed  one  at  1.219  +  0.012  X 

(40) . 

The  equilibrium  geometry  of  the  formamlde  dimer  constructed  by  Kuprlevlch  and 
Kudrltskaya  (24)  Is  given  In  the  column  headed  by  ST0-3G  +  Obs.  The  Internal 
coordinates  for  the  formamlde  molecules  In  the  dimer  were  taken  from  the 
experimental  Investigation  of  the  formamlde  monomer  by  Kurland  and  Wilson  In  1957 

(41) .  Although  there  Is  a  slight  difference  In  the  Cj=0i  bond  length,  Kurland  and 
Wilson's  (41)  results  agree  with  the  more  recent  spectroscopic  study  by  Hirota,  et. 
al  In  1974  (rightmost  column  In  Table  I)  (40).  The  hydrogen  bond  length  (rQ  ...H4) 

calculated  by  Kuprlevlch  and  Kudrltskaya  (1.853  %)  (24)  and  the  present  work  (1.651 
A)  are  essentially  Identical.  However,  we  calculate  the  0j — H4“N2  bond  an91e  t0  be 

slightly  bent  at  175°  while  Kuprlevlch  and  Kudrltskaya  (24)  constrained  It  to  be 
180  .  This  difference  should  not  be  too  significant*  because  the  potential  energy 
curve  associated  with  the  bending  of  hydrogen  bonds  Is  very  shallow  (24). 

The  Nonresonance  Excl ton- Phonon  Interaction  Parameter*  Xj* 

The  calculation  of  X^  by  Kuprlevlch  and  Kudrltskaya  (24)  and  the  present  work 
can  be  followed  by  reading  Table  II  from  left  to  right.  Note  that  R  In  Table  II  is 
equivalent  to  Rm>  ^  in  the  text.  We  see  that  the  values  of  K  calculated  by 


Kuprlevlch  and  Kudrltskaya  (34)  (1518  Nt/m  at  AR=0.0  Bohr  and  1428  Nt/m  at 
AR=0 .66653  Bohr)  are  less  than  that  of  the  present  work  (1585  Nt/m  atAR=0.0  Bohr 
and  1578  Nt/m  at  6R=0. 66653  Bohr).  Furthermore*  the  change  In  K  with  respect  to  the 
change  In  R  Is  calculated  to  be  much  greater  by  Kuprlevlch  and  Kudrltskaya  (-90 
Nt/m)  (24)  than  the  present  work  (-7  Nt/m).  In  other  words.  Kuprlevlch  and 
Kudrltskaya  (24)  have  calculated  a  greater  perturbative  effect  of  the  motion  of  the 
atom  (the  acoustical  phonon)  upon  the  Ci=0i  bond  stretching  vibration  (the  amlde- 
I  exclton)  than  Is  determined  by  our  computations.  This  result  Is  summarized  In 
Table  III  which  reveals  the  approximately  one  order  of  magnitude  difference  between 
the  nonresonance  exclton- phonon  coupling  parameter  calculated  by  Kuprlevlch  and 
Kudrltskaya  (24)  (-2.65  to  -3.64  x  10_1INt)  and  the  present  work  (-0.20  to  -0.27  x 
10“nNt). 

D.  CONCLUSIONS 

According  to  Scott  (18)  Davydov  sollton  formation  In  an  alpha-helix  Is  possible 
only  when  the  magnitude  of  the  nonresonance  exclton-phonon  Interaction  parameter 
exceeds  a  threshold  equal  to  3.5  x  10“^Nt.  The  numerical  analysis  of  Kuprlevlch 
and  Kudrltskaya  (24)  shows  that  the  magnitude  of  X.+  can  exceed  such  a  threshold 
(see  Table  III).  This  result  would  provide  support  of  the  existence  of  soli  tons  In 
polypeptide  chains.  However,  our  calculated  magnitude  of  X,+  Is  an  order  of 
magnitude  less  than  that  derived  by  Scott  (18),  and  It  would  argue  against  the 
existence  of  soli tons  In  polypeptide  chains. 

The  disagreement  between  Kuprlevlch  and  Kudrltska^a's  (24)  calculated  value  for 
X,+  and  ours  was  disturbing,  and  so  we  recalculated  Xj  using  the  exact  procedure 
of  Ref  (24).  We  found  that  Kuprlevlch  and  Kudrltskaya  were  In  error.  Instead  of  3 
x  lOljjNt,  their  calculated  magnitude  of  Xj+  should  be  0.7  x  10-11Nt.  _£]though  0.7 
x  10  Nt  Is  larger  than  our  original  calculated  magnitude  of  0.2  x  10  Nt,  It 
still  Is  less  than  the  minimum  value  of  3.5  x  10-11Nt  predicted  by  Scott  (18)  to 
sustain  a  Davydov  sollton.  Thus,  neither  the  ST0-3G  calculation  In  Ref.  (24)  nor 
ours  supports  the  existence  of  Davydov  solltons  In  hydrogen-bonded  polypeptide 
chains. 

A  more  sophisticated  ab-initio  -  SCF-M0  calculation  employing  the  4-31G  basis 
set  (42)  yielded  a  magnitude  of  1.3  x  lO'^Nt  for  Xi+  (see  Appendix  I).  This 
magnitude  Is  less  than  Scott’s  minimum  value  (18),  but  It  has  Increased  relative  to 
the  one  calculated  using  the  STO-3G  basis  set.  Therefore,  as  one  Improves  the 
sophistication  of  the  ab-initio  -SCF-M0  calculation,  one  would  expect  X  .+to 
to  Increase.  It  Is  necessary  to  Improve  the  sophistication  of  the  ab-  JP.tt.lP  -  SCF- 
M0  procedure,  because  of  the  complexity  of  the  hydrogen  bond  Involved  in  the 
exclton-phonon  Interaction  (see  Appendix  I).  As  a  result,  the  question  concerning 
the  presence  of  Davydov  solltons  In  alpha-helical  polypeptides  still  can  not  be 
answered  definitively  on  the  basis  of  the  quantum  mechanical  calculations  performed 
to  date. 
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Table  I. 
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TABLE  I.  Calculated  and  observed  equilibrium  Internal  coordinates  and  ground  state 
dipole  moments  for  the  formamlde  d1mera  and  monomer 


1.014  0.995 


TABLE  I.  (cont) 


a)  The  structure  of  the  formamlde  dimer  and  the  numbering  of  the  atoms  Is  given  In 
Figure  18b.  The  formamlde  monomer  Is  defined  by  the  molecular  fragment 
composed  of  atoms  0^»  C^,  N^»  Hj,  and  In  Figure  18b. 

b)  r  5  bond  length  between  atoms  x  and  y;  Zxy,5.  angle  between  the  x-y  and  y-z 
bonds;  l£|  =  magnitude  of  the  ground  stateya1pole  moment; 

angle  between  the  Cj-Ni  bond  vector  and  the  ground  state  dipole  moment  vector. 
Note  that  na  =  not  available. 

c)  This  work.  The  geometry  of  the  formamlde  dimer  was  determined  by  minimizing 
the  energy  of  the  molecular  system  with  respect  to  the  atomic  coordinates.  The 
energy  was  calculated  using  the  ab  Initio  self-consistent  field  (SCF)  molecular 
orbital  (MO)  method  with  the  STO-3G  minimal  atomic  orbital  basis  set.  Planar 
means  that  all  the  atoms  In  the  dimer  were  constrained  to  remain  In  one  plane 
during  the  optimization  calculation. 

d)  Ref.  24.  The  geometry  of  the  formamlde  dimer  was  determined  In  two  parts:  (1) 
the  Internal  coordinates  of  the  two  formamlde  monomer  units  were  taken  to  be 
Identical  and  are  the  observed  values  for  the  formamlde  monomer  given  in  Ref. 
41;  (2)  the  energy  of  the  dimer  was  minimized  with  respect  to  the  r(0j...H4) 
distance*  where  the  C^,  Oi*  and  N2  atoms  were  constrained  to  be  on  one 
line.  This  energy  was  calculated  using  the  ab  initio  SCF  MO  method  with  the 
ST0-3G  basis  set. 

e)  This  work.  See  c)  with  the  exception  that  dimer  is  replaced  by  monomer. 

f)  Ref.  37.  The  geometry  of  the  formamlde  monomer  was  determined  by  minimizing 
the  energy  of  the  molecule  with  respect  to  the  atomic  coordinates.  The  energy 
was  calculated  using  the  ab  Initio  SCF  MO  method  with  the  ST0-3G  basis  set. 
The  monomer  was  not  constrained  to  remain  planar  in  this  calculation,  and  It 
was  found  that  the  minimum  energy  geometry  was  not  planar  with  the  Oj-Cj-Nj-H^ 

and  01-Cj-Nj-H2  dihedral  angles  equal  to  146.8°  and  22. 8°,  respectively. 

g)  Ref.  37.  The  geometry  of  the  formamlde  monomer  was  determined  by  minimizing 
the  energy  of  the  molecule  with  respect  to  the  atomic  coordinates.  The  energy 
was  calculated  using  the  ab  Initio  SCF  MO  method  with  the  4-3 1G  spl It- valence 
basis  set.  The  monomer  was  calculated  to  be  planar. 

h)  Ref.  40.  A  microwave  spectrum  of  the  formamlde  monomer  In  the  gaseous  state. 
The  monomer  was  determined  to  be  planar. 


38 


TABLE  II.  (cont) 


a)  The  structure  of  the  formamlde  dimer  Is  given  In  Figure  18b.  The  non-resonance 
excl ton- phonon  coupling  parameter  Is  defined  In  h)  of  Table  III. 

b)  The  change  In  hydrogen  bond  length  (rlO^-fy)  In  Figure  18b)  with  respect  to  Its 
equilibrium  value.  A  positive  AR  reflects  an  Increase  In  r(Oj-H4). 

c)  The  magnitude  of  the  change  In  the  carbonyl  bond  length  (r(C.*0j)  in  Figure 
18b)  with  respect  to  Its  equilibrium  value.  Note  that  a  post tl ve/negative  f 
reflects  an  Increase/decrease  of  rtC^Oj). 

d)  The  total  ground  electronic  state  energy  of  the  dimer. 

e)  The  present  work  Is  represented  by  P.  The  ab  Initio  self-consistent  field 
(SCF)  molecular  orbital  (MO)  method  with  the  ST0-3G  minimal  atomic-orbital 
basis  set  was  used  to  calculate  the  total  energy  of  the  dimer  as  a  function  of 
AR  and  Jf  .  Note  the  following:  (1)  J  was  positive*  as  It  was  increased  In 
Increments  of  0.02;  12)  5  was  Increased  by  moving  the  Oj  atom  and  the  molecular 
fragment  composed  of  the  C^.  N^,  Hj»  H2*  and  atoms  (see  text). 

f)  The  work  of  Kuprlevlch  and  Kudrltskaya  In  Ref.  24  Is  represented  by  KK.  The  At 

Initio  SCF  MO  method  with  the  ST0-3G  basis  set  was  used  to  calculate  the  total 

energy  of  the  dimer  as  a  function  of  &R  and  5  .  Note  the  following:  (1)  $ 

was  negative*  as  It  was  decreased  In  Increments  of  0.02;  (2)  $  was  decreased  by 
moving  only  the  0j  atom  (see  text);  (3)  the  total  energy  of  the  dimer  decreases 
as  J  decreases,  because  the  Initial  r(C^*0i)  *  1.243  X  Is  greater  than  the 
optimal  r(Ci=0i)  *  1.219  X  calculated  using  the  ab-initio  SCF  MO  method  with 
the  ST0-3G  basis  set  (see  text  and  Table  I). 

g)  The  change  In  total  energy  associated  with  the  change  In  5  (A?) 

h)  The  change  In  AE  associated  with  Af  . 

I)  The  estimate  of  the  diagonal  force  constant  for  the  stretching  of  the  C^O^ 

bond.  K  2s  A( A  E ) /  (A?  )2. 

J)  The  third  derivative  of  the  total  energy  of  the  dimer  which  Is  the  crucial 
component  of  the  non-resonance  excl ton- phonon  coupling  parameter  (see  text  and 
Table  III). 
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TABLE  III.  Calculated  values  of  the  non-resonance  excl ton- phonon 

coupling  parameter  for  different  diagonal  force  constants  for 
the  C  *  0  bond  In  the  formamlde  dimer®  and  monomer8 1 


1150  -0.27  -3.64 

(monomer.  Ref.  40)c^ 


1380  -0.22  -3.03 

(monomer,  Ref.  37)^ 

1540  -0.20  -2.72 

(monomer.  Ref.  24)®^ 

1473  -0.21  -2.84 

(dimer.  Ref.  24)f) 

1581 

(dimer.  This  work)95 


0.20 


2.65 


TABLE  III.  (cont) 


a)  The  structure  of  the  formamlde  dimer  Is  given  In  Figure  18b.  The  formamlde 
monomer  Is  defined  by  the  molecular  fragment  composed  of  atoms  Op  Cp  Np  H 
H2,  and  H3  in  Figure  18b. 

b)  The  diagonal  force  constant  for  the  stretching  of  the  bond. 

c)  Derived  from  spectroscopic  data. 

d)  Calculated  using  the  ab-initio  self-consistent  field  (SCF)  molecular  orbital 
(MO)  method  with  the  4-31G  split-valence  basis  set. 

e) -g)  Calculated  using  the  ab-initio  SCF  MO  method  with  the  ST0-3G  basis  set. 
h)  The  non-resonance  excl ton- phonon  coupling  parameter. 

Xl+”  ^  '5^'  *here  *  the  experimental  mean  frequency  of  the 

amide- 1  vibration  (1660  cm”*  *  3.297  x  10”^J)» 

K  as  defined  in  b)«  and  AK/AR  as  defined  In  1)  of  Table  II. 

I)  The  value  of  Xj*  determined  using  &K/AR  calculated  In  the  present  work. 

J)  The  value  of  X.*  determined  using  4K/AR  calculated  by  Kuprlevlch  and 
Kudrltskaya  In  Ref.  24. 
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APPENDIX  I.  The  calculation  of  the  nonlinear  coupling  constant  (x/),  using  an 
ab-initio  -SCF-MO  procedure  with  the  4-316  basis  set/ 

The  calculation  of  X^+  using  an  ab-initio  -SCF-MO  formalism  with  the  4-31G 
split  valence  basis  set  Is  a  great  Improvement  over  earlier  calculations  employing 
the  ST0-3G  basis  set  (see  Section  III  and  Ref.  1-1).  However*  before  discussing  the 
results  of  the  4-31G  calculation,  there  Is  one  finding  to  be  noted  concerning  the 
ST0-3G  calculations;  the  magnitude  of  «  3  x  10-11Nt  calculated  by  Kuprlevlch  and 

Kudrltskaya  ( I— 1 )  Is  too  high.  We  recalculated  X,+  by  following  exactly  the 
procedure  outlined  by  Kuprlevlch  and  Kudrltskaya  tl-1),  and  obtained  a  magnitude 
equal  to  0.7  x  KT^Nt.  Although  Kuprlevlch  and  Kudrltskaya  (1-1)  mention  that  they 
had  SCF  convergence  problems  (we  had  none),  the  errors  In  their  calculation  of  X^+ 
are  not  clear  to  us.  In  conclusion,  neither  the  ST0-3G  calculation  In  Ref.  1-1  nor 
ours  (Section  III)  yields  a  value  of  \X,+|wh1ch  exceeds  the  threshold  of  3.5  x  10"11 
Nt  predicted  by  Scott  (1-2).  This  result  will  be  discussed  In  detail  below. 

Ab-Initio  -SCF-MO  calculations  of  molecular  properties  employing  the  ST0-3G 
basis  set  are  crude.  This  fact  Is  particularly  true  when  treating  Interactions  like 
Xx+,  which  are  associated  with  hydrogen  bonding.  The  starting  point  for  the 
treatment  of  hydrogen-bonded  systems  Is  to  use  the  4-31G  basis  set.  It  Is  stressed 
that  this  basis  set  Is  only  the  starting  point  because  the  van  der  Waals  forces 
which  are  important  in  hydrogen  bonding  are  not  being  described  reasonably  until  one 
utilizes  second  and  third  order  M^ller-Plesset  perturbation  theory  In  conjunction 
with  polarized  basis  sets  of  the  6-311G**  type  (1-3). 

The  ab-initio  -SCF-MO  calculation  of  X1+  using  the  4-31G  basis  set  followed  the 
same  procedure  outlined  for  the  calculation  employing  the  ST0-3G  basis  set:  (1)  use 
a  formamlde  dimer  as  the  simplest  model  of  the  polypeptide  chain,  (2)  calculate  an 
equilibrium  geometry  of  the  dimer,  (3)  calculate  K  or  K^q  as  a  function  of 
or  ro...H*  Ho*ever,  unlike  the  STQ-3G  calculation  which  determined  Kc=0  at  two 
values'of  Rq  *  0.00  and  0.35  A,  the  4-31G  calculation  computed  Kq_q  at  seven 

values  of  Rq. ..H  *  -  °‘05,  ±  0.02,  ±  0.01,  and  0.00  X. 

Calculated  (STO-3G  and  4-3 1G)  and  observed  equilibrium  geometries  of  the 
formamlde  monomer  are  shown  In  the  top  part  of  Figure  1-1.  Note  that  the  calculated 
geometries  were  constrained  to  be  planar,  and  the  observed  one  was  found  to  be 
planar.  It  Is  clear  that  the  4-31G  calculated  geometry  matches  up  better  with  the 
observed  one  than  does  the  ST0-3G  calculated  geometry.  This  observation  Is 
particularly  true  for  the  Cj-N^  bond  length.  The  term  |£|  Is  the  magnitude  of  the 
dipole  moment,  and  6 {ft,  Jc^-Nj^  5  0  Is  the  angle  between  the  dipole  moment  vector 
and  Cj-N^  bond  vector. 

Calculated  (ST0-3G  and  4-3 1G)  equilibrium  geometries  of  the  formamlde  dimer  are 
shown  In  the  bottom  half  of  Figure  1-1.  No  observed  equilibrium  geometry  of  the 
dimer  Is  available  In  this  configuration.  Similar  to  that  for  the  monomer,  the 
calculated  geometry  of  the  dimer  was  constrained  to  be  planar.  The  major  points  to 
be  noted  concerning  the  calculated  dimer  geometries  are  (1)  the  Rq  *0}  and  RH  -n2  In 

the  dimer  with  respect  to  Rq  bq  and  RH  _N  In  the  monomer  are  greater  for  the  4-3 1G 

11q  11 

calculation  (+  0.005  and  +  0.007  l,  respectively)  than  those  for  the  ST0-3G 
calculation  (♦  0.001  and  +  0.003  X,  respectively).  Clearly,  the  formamlde  monomer 
Is  calculated  In  the  4-3 1G  basis  set  to  be  more  affected  by  hydrogen  bonding  than  in 
the  STO-3G  basis  set. 
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.  Equilibrium  geometries  of  the  formamlde  monomer  and  dimer 

calculated  using  the  ab-initio  -SCF-MO  procedure  with  the  ST0-3G 
and  4-31G  basis  sets.  The  observed  geometry  of  the  formamlde 
monomer  determined  by  HI  rota*  fit  .  a1  (1-6)  Is  also  shown. 
Note  that  \/t  I  Is  the  magnitude  of  the  molecular  dlple  moment* 
and  ©  (yw,»  rc  _N  >  Is  the  angle  between  the  dipole  moment 

vector  and  the  Cj-Nj  bond  vector. 
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Calculated  (ST0-3G  and  4-31G)  electronic  charge  distributions  for  the 
calculated  equilibrium  geometries  of  the  formamlde  monomer  are  shown  In  the  left- 
half  of  Figure  1-2.  The  polarization  of  charges  across  all  bonds  Increases  as  one 
goes  from  the  ST0-3G  to  the  4-3 1G  basis  set.  It  Is  also  shown  that  the  charge 
distribution  for  the  formamlde  monomer  calculated  using  the  even  more  elaborate  6- 
31G  basis  set  (1-4)  Is  not  too  different  from  the  4-31G  calculated  charge 
distribution.  This  statement  Is  supported  further  by  noting  that  Ifc,  I  (4-31G)  = 
4.471D  Is  closer  In  value  to  I H  l  (observed)  =  3.71  +  0.06D  than  is  l/£  1(ST0-3G)  = 
2.641D.  Both  6  (4-31G)  *  41.3°  and  0  (ST0-3G)  *  39.3°  are  In  near-agreement  with  6 
(observed)  *  39.6°. 

Calculated  (STO-3G  and  4-31G)  electronic  charge  distributions  for  the 
calculated  equilibrium  geometries  of  the  formamlde  dimer  are  shown  In  the  right-half 
of  Figure  1-2.  As  In  the  case  of  the  formamlde  monomer,  the  polarization  of  charges 
across  all  bonds  Increases  as  one  goes  from  the  ST0-3G  to  the  4-31G  basis  set.  This 
observation  Is  reflected  In  f/5 (  (4-31G)  =  10.331D  being  greater  than  |  ($ TO- 
3G)  *  6.280D.  The  greater  polarization  of  charge  across  the  hydrogen  bond 
calculated  using  the  4-31G  basis  set  Implies  a  greater  Intermolecular  Interaction 
which  should  result  In  a  greater  value  of  |Xj\ 

Figure  1-3  presents  the  plot  of  Kq_q  Rq  >h  calculated  with  the  4-31G  basis 
set.  The  linear  regression  fit  of  the  points  yielded  an  equation,  Kq_0  =  1463.1  - 

97.5  Rq  with  an  r^  *  0.993.  The  slope  of  the  plot  yields  a  value  of  Xj  +  = 

-0.13  x  10~10Nt.  The  terms^Xj*,  (force  constant  for  the  calculated  equilibrium 

geometry  at  Rq  H  *  1.922  A),  tiw  iamlde-1  frequency),  ESCF  (total  energy  of  the 
dimer),  and  A  I61nd1ng  energy  of  the  dimer  or  hydrogen  bond  stabilization  energy) 
calculated  for  the  formamlde  dimer  employing  the  4-31G  and  ST0-3G  basis  sets  are 
given  In  Table  I-I. 

Although  Ix^l  (4-31G)  has  Increased  relative  to  |Xj+|  (ST0-3G),  It  Is  still 

less  than  the  minimum  value  of  0.35  10"*®Nt  predicted  by  Scott  (1-2)  to  sustain  a 

Davydov  sol  1 ton  In  alpha-helical  polypeptides.  However,  since  |X^+|  Is  proportional 
to  the  hydrogen  bond  binding  energy  (Ain  Table  I-I),  the  stronger  the  hydrogen  bond, 
the  greater  the  value  of  lx^+|  .  There  Is  no  doubt  that  the  hydrogen  bond  for  the 
formamlde  dimer  is  stronger  than  what  we  have  calculated  using  the  4-31G  basis  set, 
because  we  have  neglected  the  van  der  Weals  forces.  Furthermore,  the  strength  of 
hydrogen  bonds  between  molecular  units  In  a  hydrogen-bonded  polymer  chain  can  be 
double  that  of  the  hydrogen  bond  In  a  dimer  (1-5).  Therefore,  by  Including  van  der 
Waals  forces  In  our  theoretical  formalism  and  by  Increasing  the  number  of  formamlde 
monomer  units  In  our  model  of  the  polypeptide,  we  would  expect  the  calculated  |Xj*| 

to  increase  to  a  value  perhaps  greater  than  0.35  x  10”*°Nt.  From  the  perspective  of 
the  quantum  mechanical  calculation  of  Xj  ,  the  question  of  the  presence  of  Davydov 
solltons  in  alpha-helical  polypeptides  Is  still  not  resolved. 


Figure  1-2.  Electronic  charge  distributions  for  equilibrium  geometries 
of  the  formamlde  monomer  and  dimer  calculated  using  the  ab- 
Inltlo  -SCF-MO  procedure  with  the  ST0-3G  and  4-31G  basis  sets. 
The  electronic  charge  distribution  for  the  formamlde  dimer 
calculated  by  Hagler  and  Laplcclrella  (1-4)  with  the  6-31G  basis 
set  Is  also  shown.  The  terms  |/M{  and  0  i£,  Jq  _n  )  are 

defined  In  the  caption  to  Figure  1-1.  11 
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TABLE  1-1 


Calculated 
parameters  for 

values  of  the  exclton-phonon  coupling  term  (X^+)  and 
the  formamlde  dimer. 

PARAMETER 

PRESENT  WORK 

4-3 1G 

PRESENT  WORK 
ST0-3G 

SCOTT,  1982 

*  x+  ao10Nt) 

-0.13 

=£LQ2 

AJL6 

Kc?^(Nt/m) 

1275c) 

1553c) 

na 

« 

n 

3 

1 

M 

1660 

1660 

na 

ESCF(au) 

-337.376457664 

-333.383404267 

na 

A  *le\) 

-0.36168 

-0.18984 

na 

(a)  A  *  £$CF  (dimer)  -2ESCp  (monomer);  na  *  not  applicable 

(b)  A.C.  Scott,  Phys.  Rev.  A  26  ,  578  (1982) 

(c)  Kc?£  (4-31G.  monomer)  *  1285  Nt/m; 

K  ®9(STO-3G»  monomer)  *  1586  Nt/m. 
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